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ABSTRACT 

We present a series of cosmological TV-body simulations which make use of the hy- 
drodynamic approach to the evolution of structures (Dommguez 2000). This approach 
addresses explicitly the existence of a finite spatial resolution and the dynamical effect 
of subresolution degrees of freedom. We adapt this method to cosmological simulations 
of the standard ACDM structure formation scenario and study the effects induced at 
redshift z = by this novel approach on the large-scale clustering patterns as well as 
(individual) dark matter halos. 

Comparing these simulations to usual A/-body simulations, we find that (i) the 
new (hydrodynamic) model entails a proliferation of low-mass halos, and (ii) dark 
matter halos have a higher degree of rotational support. These results agree with 
the theoretical expectation about the qualitative behaviour of the " correction terms" 
introduced by the hydrodynamic approach: these terms act as a drain of inflow kinetic 
energy and a source of vorticity by the small-scale tidal torques and shear stresses. 

Key words: gravitation - methods: numerical - methods: TV-body simulations - 
galaxies: formation - cosmology: theory 



1 INTRODUCTION 

Gravitational instability is commonly accepted as the basic 
mechanism for structure formation on large scales. Com- 
bined with the CDM model it leads to the picture of hier- 
archical clustering with wide support from deep galaxy and 
cluster observations. During the recent phase of cosmic evo- 
lution groups and clusters of galaxies condense from large 
scale density enhancements, and they grow by accretion and 
merger processes of the environmental cosmic matter. 

But despite the fact that the currently favoured ACDM 
model has proven to be remarkably successful on large scales 
(cf. WMAP results, Spergel et al. 2003, Spergel et al. 2006), 
recent high-resolution A-body simulations still seem to be in 
contradiction with observation on sub-galactic scales. There 
is, for instance, the problem with the steep central densities 
of galactic halos as the highest resolution simulations fa- 
vor a cusp with a logarithmic inner slope for the density 
profile of approximately —1.2 (Diemand, Moore & Stadel 
2005; Fukushige, Kawai & Makino 2004; Power et al. 2003), 
whereas high resolution observations of low surface bright- 
ness galaxies are best fit by halos with a core of constant 
density (Simon et al. 2005; de Block & Bosma 2002; Swa- 
ters et al. 2003). A further problem relates to the overabun- 



dance of small-sized (satellite) halos; there are many more 
subhaloes predicted by cosmological simulations than actu- 
ally observed in nearby galaxies (e.g., Moore et al. 1998, 
Klypin et al. 1999, Gottlober et al. 2003). The lack of obser- 
vational evidence for these satellites has led to the suggestion 
that they are completely (or almost completely) dark, with 
strongly suppressed star formation due to the removal of gas 
from the small protogalaxies by the ionising radiation from 
the first stars and quasars (Bullock et al. 2000; Tully et al. 
2002; Somerville 2002; Hoeft et al. 2005). Other authors sug- 
gest that perhaps low mass satellites never formed in the 
predicted numbers in the first place, indicating problems 
with the ACDM model in general, which is replaced with 
Warm Dark Matter instead (Knebe et al. 2002; Bode, Os- 
triker & Turok 2001; Colin et al. 2000). Suggested solutions 
also include the introduction of self-interactions into colli- 
sionless A-body simulations (e.g. Spergel & Steinhardt 2000; 
Bento et al. 2000), and non-standard modifications to an 
otherwise unperturbed CDM power spectrum (e.g. bumpy 
power spectra: Little, Knebe & Islam 2003; tilted CDM: Bul- 
lock 2001c). Recent results from (strong) lensing statistics 
though suggest that the predicted excess of substructure is 
in fact required to reconcile some observations with theory 
(Dahle et al. 2003, Dalai & Kochanek 2002), but this con- 
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elusion has not been universally accepted (Sand et al. 2004; 
Schechter & Wambsganss 2002; Evans & Witt 2003). 

The discovery of the mismatch between observations 
and simulations is a result of the increase in the resolu- 
tion of TV-body simulations over the last years. This has 
emphasized the importance of the influence of subresolution 
scales on the simulated dynamics, at least when it comes 
to halo properties. The purpose of this work is to study 
the hydrodynamic-like formulation of the formation of cos- 
mologial structure proposed recently by Dommguez (2000), 
dubbed SSE (small-size expansion). The SSE addresses ex- 
plicitly the existence of a finite spatial resolution and the 
dynamical effect of subresolution degrees of freedom. Al- 
though developed independently, the SSE approach is close 
in spirit to the Large-Eddy Simulation of turbulent flow 
(see, e.g., Pope 2000 and refs. therein). This is a method 
devised to simulate only the relevant, large-scale degrees of 
freedom according to the Navier-Stokes equations describing 
flow in the high-Reynolds number (i.e., turbulent) regime: 
physically meaningful approximations are made in order to 
model the coupling to the neglected, small-scale degrees of 
freedom. 

The SSE starts from the microscopic equations of mo- 
tion for a set of TV particles under their mutual gravity 
and provides a set of hydrodynamic-like equations for the 
(coarse-grained) mass density and velocity fields. These new 
equations now contain "correction terms" which describe 
the effects of the coarse-graining procedure and correct for 
them, respectively. It therefore only appears natural to im- 
plement these correction terms into an (adaptive) mesh TV- 
body code where the density is treated in a coarse-grained 
fashion, too: mesh-based Poisson solvers frequently used for 
cosmological simulations smooth the particle distribution 
onto a grid and hence deal with a coarse-grained density 
field when solving for the potential via Poisson's equation. 
For this purpose we will adapt the open source TV-body code 
MLAPM (Multi-Level Adaptive Particle-Mesh)*. The particles 
of the TV-body simulations presented in this study are effec- 
tively hydrodynamical Lagrangian particles which move un- 
der the action not only of the mesh-computed gravitational 
force, but also of the additional, correction terms modeling 
sub-grid degrees of freedom in the context of the SSE. 

The rest of the work is structured as follows: in Sec. 2 we 
describe the theoretical background of the SSE and provide 
the link to mesh-based TV-body codes. In Sec. 3 we present 
the simulation of several models (standard ACDM model, 
AWDM model, and two runs incorporating the SSE correc- 
tions). In Sec. 4 we perform a comparative analysis of the 
four runs at redshift z — from two complementary points 
of view: properties of the mass density and velocity fields, 
and properties of halos. Finally, Sec. 5 includes a discussion 
of the results and the conclusions. 



2 THE HYDRODYNAMIC APPROACH 

We deal with a system of nonrelativistic, identical point par- 
ticles which (i) are supposed to interact with each other via 
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gravity only, (ii) look homogeneously distributed on suffi- 
ciently large scales, so that the evolution at such scales cor- 
responds to an expanding Friedmann-Lemaitre cosmological 
background, and (iii) deviations to homogeneity are relevant 
only on scales small enough that a Newtonian approxima- 
tion is valid to follow their evolution. Let a(t) then denote 
the expansion factor of the Friedmann-Lemaitre cosmologi- 
cal background, H(t) = a /a the associated Hubble function, 
and Qb(t) the homogeneous (background) density on large 
scales. x a is the comoving spatial coordinate of the a-th 
particle, u a — ax a its peculiar velocity, and m its mass. In 
terms of these variables the evolution is described by the fol- 
lowing set of equations (Peebles 1980) (V Q denotes a partial 
derivative with respect to x a ): 
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where w a is the peculiar gravitational acceleration acting 
on the a-th particle. Finally, Eqs. (1) must be subjected to 
periodic boundary conditions in order to yield a Newtonian 
description consistent with the Friedmann-Lemaitre solution 
on large scales (Buchert & Ehlers 1997). 

If we now assume that the actual measure of the density 
field in an TV-body code depends on a smoothing window 
W(z), the microscopic field g m ic relates to the measured 
(coarse-grained) field g in the following way: 



Qm ic {*,t) = *W(X-Xa(t)), (2a) 

a 

e (x,t;L) = J% 



li' ( ) g mic (y,t). (2b) 



The physical interpretation of the field g(x; L) follows 
straightforwardly from the properties of the smoothing win- 
dow: it is proportional to the number of particles contained 
within the coarsening cell of size w L centered at x. A mi- 
croscopic peculiar-momentum density field and the corre- 
sponding coarse-grained field can be defined in the same 
way: 

j mic (x,t) = ^3 5^u a (t) <5 (3) (x-x a (*:)), (3a) 

a 

j(x,t;L) = w(^^jj mic (y,t). (3b) 



One can introduce peculiar velocity fields u mic and u by 
definition as j =: gu and similarly for n m ic- The physical 
meaning of u(x; L) is also simple: it is the center-of-mass 
peculiar velocity of the subsystem defined by the particles 
inside the coarsening cell of size w L centered at x. Notice 
that u is not obtained by coarse-graining u m i C : from a dy- 
namical point of view, it is more natural to coarse-grain the 
momentum rather than the velocity, since the former is an 
additive quantity for a system of particles. Finally, one can 
define peculiar gravitational acceleration fields w m i C and w 
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through coarse-graining of the force: 

£WW mic (x,t) = ^3^w a (t) <5 (3) (x - x„(()), (4a) 



gw(x,t;L) 



(ty 

L3 l L 



w mic (y,$±>) 



The field w(x) has the physical meaning of the center-of- 
mass peculiar gravitational acceleration of the subsystem 
defined by the coarsening cell at x. 

From these definitions and Eqs. (f), it is straightfor- 
ward to derive the evolution equations obeyed by the coarse- 
grained fields g and u (from now on, d/dt is taken at con- 
stant x and L, and V means partial derivative with respect 
to x): 

^+3ffe = -V(eu), (5a) 

^4^- +4Hgu= gw- -V- (ouu + n), (5b) 
at a 

where a new second-rank tensor field has been defined 
(dyadic notation): 



n(x,t;L) = 



w |x ~ y| 

L3 I L 



Qmic{y,t) 



(6) 



[u mic (y,t) - u(x,t;L)][u mic (y,t) - u(x,t;L)]. 

The field n(x) is due to the velocity dispersion, i.e., to the 
fact that the particles in the coarsening cell have in gen- 
eral a velocity different from that of the center of mass. The 
physical meaning of Eqs. (5) is simple: they are just bal- 
ance equations, stating mass conservation and momentum 
conservation, respectively. The term g w codifies the gravita- 
tional interaction between the coarsening cells and does not 
satify, in general, Poisson's equation or the curl-free condi- 
tion. The term V • II represents a kinetic pressure force due 
to the exchange of particles between neighboring coarsening 
cells (just like for an ideal gas) and it has the same physical 
origin as the convective term V • (puu), i.e., a nonlinear 
mode-mode coupling of the velocity field. The difference is 
that the convective term couples only modes on scales > L, 
while the velocity dispersion term codifies the dynamical ef- 
fect of the coupling of the modes on scales > L with the 
modes on scales < L. Eqs. (5) are exact: as one changes the 
smoothing length, the fields g, u, w, II change in such a 
way that the form of the equations remains the same (for 
example, upon increase of the smoothing length, part of the 
dynamical effect described by V ■ (f>uu) is shifted towards 
V • II). This property is reflected in that the equations are 
not an autonomous system for g and u. In fact, they are just 
the first ones of an infinite hierarchy, as can be checked by 
computing the evolution equations for the fields w and II. 
To obtain a useful set of equations, it is necesary to truncate 
this hierarchy by looking for a functional dependence of w 
and II on g and u. 

The peculiarities of the problem at hand (collisionless 
matter in the non-stationary state of structure formation) 
prevent the usual truncation of the hierarchy leading to 
the Euler or Navier-Stokes equations, respectively (see, e.g., 
Chapman & Cowling 1991). The small-size expansion (SSE) 
is a specific truncation for this problem (Dommguez 2000, 
2002; Buchert & Dommguez 2005), that starts from the 
physical assumption that the coupling to the small-scales 



is weak (this can be argued on the basis that, in a hierarchi- 
cal scenario, the smaller scales "virialize" earlier and thus 
"decouple" from the evolution of the larger scales). Then 
the fields II and w are derived as a formal expansion in L: 
Keeping terms up to order L 2 , Eqs. (5) become {di = d/dxi\ 
summation over the repeated index i is understood) 
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with the additional acceleration 
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The constant B is determined by the smoothing window 

W(z), 



dz z 2 W(z) : 
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dz z 4 W(z). 
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To order L°, Eqs. (7) reduce to the "dust (pressureless) ap- 
proximation" for cosmological structure formation (Sahni 
& Coles 1995): w m represents the mean-field gravity cre- 
ated by the monopole moment of the matter distribution 
in the coarsening cells, i.e., the total mass, and the cou- 
pling to the small-scales is neglected altogether. To order 
L 2 there are two kinds of corrections: the term proportional 
to the tidal tensor Vw mf (stemming from a term w — w mf in 
Eqs. (5)) models the gravitational force of the higher-order 
multipole moments, i.e., the coupling to the subresolution 
configurational degrees of freedom; the term proportional to 
(9ju)(9ju) (stemming from II in Eqs. (5)) models the effect 
of velocity dispersion, that is, the coupling to the subreso- 
lution kinetic degrees of freedom. 

The expected dynamical effect of these new terms with 
respect to the "dust evolution" has been studied theoreti- 
cally (Dommguez 2000, 2002; Buchert & Dommguez 2005). 
There is evidence that, assuming a locally plane-parallel 
collapse, these terms mimick the "adhesion model" (Gur- 
batov, Saichev & Shandarin 1989, Sahni & Coles 1995), in 
which recently collapsed regions stabilize — or more gen- 
erally speaking, the term due to velocity dispersion tends 
to reduce the inflow velocity in collapsing regions and fa- 
vors the formation of gravitationally bound systems. It has 
also been shown that the correction terms act as a source 
of vorticity via small-scale tidal forces and shear stresses. 
The "dust model" lacks a source of vorticity w = V x u, 
and the initially present one is damped by the cosmological 
expansion in the linear regime. Thus, in that respect the cor- 
rections to the "dust model" can be particularly important. 
Taking the curl of Eq. (7b) we obtain 



9oj 1 „ . N _ 

— = - Hu + -V x (u x w) + V x C , 

at a 



(10) 



where the term V x C is a source of vorticity, i.e., it does not 
vanish in general even if us = 0, as has been confirmed per- 
turbatively by Dommguez (2002). We comment further on 
the relationship between vorticity and angular momentum 



© 0000 RAS, MNRAS 000, 000-000 



4 Knebe A. et al. 



in App. A, where also some results concerning the conser- 
vation of energy, momentum, and angular momentum are 
derived. 

The dynamical evolution predicted by Eqs. (7) can 
be implemented without much difficulties in a particle- 
mesh (PM) code of TV-body simulation. Mass conservation, 
Eq. (7a), is automatically satisfied by the code. The acceler- 
ation w mf given by Eqs. (7c) and (7d) agrees with the value 
returned by the Poisson solver on a grid, and the grid con- 
stant sets naturally the resolution L. In principle, one only 
needs to take care of Eq. (7b), which can be re- written in 
Lagrangian coordinates as 

u = w mf -ffu + C. (If) 

The computation of C, given by Eq. (8), on the grid of the 
Poisson solver is a highly non-trivial but manageable task. 
Thus, Eq. (11) — together with x = u/a, see Eq. (la) - 
determines the motion of Lagrangian fluid elements, which 
are sampled by the particles of the simulation. If we set 
B = (=> C = 0) in Eq. (11) we recover the equations of 
motions that are being integrated in a standard PM code 
for the update of particle velocities and positions during the 
course of the simulation. 



3 THE iV-BODY SIMULATIONS 
3.1 The Setup 

The JV-body simulations presented in this study were carried 
out using a version of the open source adaptive mesh refine- 
ment code MLAPM (Knebe, Green & Binney 2001). This code 
reaches high force resolution by refining high-density regions 
with an automated refinement algorithm. These adaptive 
meshes are recursive: refined regions can also be refined, 
each subsequent refinement having cells that are half the 
size of the cells in the previous level. This creates a hier- 
archy of refinement meshes of different resolutions covering 
regions of interest. The refinement is done cell by cell (in- 
dividual cells can be refined or de-refined) and meshes are 
not constrained to have a rectangular (or any other) shape. 
The criterion for (de-) refining a cell is simply the number of 
particles within that cell and a detailed study of the appro- 
priate choice for this number can be found elsewhere (Knebe 
ct al. 2001). The code also uses multiple time steps on dif- 
ferent refinement levels where the time step for each level 
is two times smaller than the step on the previous level. 
The latest version of MLAPM also includes an adaptive time 
stepping that adjusts the actual time step after every major 
step to restrict particle movement across a cell to a partic- 
ular fraction of the cell spacing, hence, fine tuning accuracy 
and computational time. 

As outlined above, the only necessary modification re- 
quired to model the "Hydrodynamic Approximation" (or 
"HAPPI" here afterwards) is to account for the correction 
term C in Eq. (11) when updating the particle velocities^. 
MLAPM has therefore been modified to not only calculate 
the density field on its hierarchy of nested refinement grids 
but also the velocity field. The V-operator and spatial 

t The modifications are part of MLAPM vl.4 (and all later versions) 
and can be switched on using -DHAPPI upon compile time. 



derivatives, respectively, have been approximated by finite- 
differences using the two nearest neighbors (in each dimen- 
sion) to the cell for which the correction term is being cal- 
culated. Cells close to a refinement boundary for which not 
enough surrounding nodes are present, obtain their correc- 
tion values interpolated downwards from the next coarser 
grid level. The assignment of mass and momentum on the 
grid is done with a triangular-shaped-cloud window (Hock- 
ney & Eastwood 1988), 

r |-2 2 for |*| < 0.5, 

W ( z )=\ Ml-N) 2 for 0.5 < 1*1 < 1.5, (12) 
I otherwise, 

for which B = 1/4 according to Eq. (9). We remark that, 
due to the dynamical (de)refinement procedure of the MLAPM 
code, the resolution is space-dependent in a discrete man- 
ner, while Eqs. (7) were derived under the assumption of 
a spatially homogeneous length L. The SSE can be gener- 
alized to the case of a smoothly varying L(x) (Dommguez, 
unpublished; Dommguez 2002 contains the generalization to 
a time-dependent L), but we have neglected this additional 
complication because the fraction of particles which are in 
regions where L jumps is less than 1% during the run. Fi- 
nally, we also note that this numerical method of integrating 
the hydrodynamic equations is different from, albeit similar 
to, the Smoothed Particle Hydrodynamics method (Gingold 
& Monaghan 1977; Lucy 1977) frequently used in cosmolog- 
ical simulations involving baryonic matter. 

We ran four CDM simulations with cosmological param- 
eters in agreement with the so-called concordance model, i.e. 
n = 1/3, A,, = 2/3, a s = 0.88, h = 2/3: 

• one standard ACDM model, 

• one AWDM model (m W DM = 0.5keV), 

• two ACDM models with B=l/4 and B=l, respectively. 

Even though B is actually determined by the smoothing 
window, we also considered a four times larger value as if 
it were a free parameter of the model. This model is to be 
understood as an "academic toy model" where we hope to 
gain better insight into the effects of the correction term on 
cosmological structure formation. 

All simulations consist of N = 128 3 particles in a box of 
side length 25/i _1 Mpc (the mass of a simulation particle is 
TxlO'r'Mo), and they were started at redshift 
z = 35. The two ACDM models with B /0 are also dubbed 
ACDMhappil (B=l/4) and ACDMhappi2 (5=1). We chose 
to also run a AWDM model to allow for a more complete 
comparision of the new HAPPI models to other, alternative 
cosmologies. A more elaborate study of the AWDM model 
and Warm Dark Matter can be found in Knebe et al. (2002). 

The force resolution in MLAPM is determined by the finest 
refinement level reached throughout the run. While all four 
models applied exactly the same refinement criterion (six 
particles per cell), the ACDMhappi2 run only invoked five 
refinement levels whereas all other runs used seven levels 
at redshift z = 0. In terms of force resolution this trans- 
lates to 10fe -1 kpc spatial resolution for ACDMhappi2 (cor- 
responding to an estimated maximum density ~ 5 x 10 4 Qb) 
and 2.5h~ kpc for all the other models (maximum den- 
sity ~ 3 x 10 6 f?(,)- This difference can be ascribed to the 
smoothening effect of the correction terms in the dynamical 
equations that are more effective for higher values of B. 
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ACDM 

ACDM happi 1 
ACDM happi 2 
AWDM 




Table 1. Fraction of HAPPI particles in the subset of low-density 
particles, and in the subset of high-density particles, respectively. 
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Figure 1. Variation with the expansion factor a of the dimen- 
sionless quantity X defined by Eq. (14). 



3.2 Accuracy of the code 

A useful check of the accuracy of a simulation is provided 
by the global invariants of the dynamical system, which we 
derive and discuss in App. A. 

The dynamical evolution conserves the quantity 
a u a (if the force and density interpolation schemes are 
identical). It was found that departures from the initially 
vanishing value satisfy the bound 



1 

JV 



< 3x lO"' 



N 

- 1 Vi 
N t-^ 

a = l 



,(13) 



where at redshift z = the average particle velocity was 
"Umean ~ 300 km s -1 in all four models. 

In standard TV-body codes it is common practice to 
check constancy of the invariant I (see Eq. (A5)) following 
from the Layzer-Irvine equation (e.g., Knebe et al. 2001). 
This is not the case for the two HAPPI models, since the 
correction term C is a source or drain of energy as discussed 
in the App. A. We though chose to plot in Fig. 1 the dimen- 
sionless quantity 



(14) 



as a function of cosmic expansion factor a, where (7 m is 
the mean-field potential energy defined in Eqs. (Al). For 
the ACDM and AWDM models, the Layzer-Irvine equation 
holds and predicts! = 0. Departures from this result are as- 
cribed to both integration/truncation errors introduced by 
the code and the fact that the particle shape is neither con- 
stant in time nor space due to the adaptive nature of MLAPM 
(cf. Knebe et al. 2001). ACDMhappil performs rather simi- 
lar to the standard CDM model, indicating that the effect of 
the correction term C in the evolution of X is small compared 
to the numerical errors. ACDMhappi2, on the other hand, 
departs noticeably from the other models and T changes sign 
at around a redshift of z ~ 3, which lets us expect to find 
larger differences between ACDMhappi2 and ACDM. 



redshift 


low-density 


high-density 
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14% 
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3.3 The Importance of the HAPPI correction 

In order to test the importance of the correction term Eq. (8) 
we calculated the ratio of the mean-field acceleration (i.e. 
F — |w mf j) and the additional acceleration (C = |C|) for 
each individual particle as a function of the local density at 
various redshifts. The result for the ACDMhappi2 model, for 
which the effect of C is the largest, can be viewed in Fig. 2. 
This figure indicates that the effect of the newly added terms 
is generally rather small especially at late times. At redshift 
z = the fraction of all particles with F/C < 1 ("HAPPI 
particles") is a mere 9% while it increases to 15% at z — 4. In 
order to examine a possible trend with density, we consider 
at each redshift two subsets of particles according to whether 
the density is above or below the virial overdensity A v i r (cf. 
Eq. (17) in Section 4.2) and hence particles of the high- 
density subset either already belong to virialized structures 
or will be part of them at a later time. Table 1 gives the 
fraction of HAPPI particles in each subset. 

This raises the question about the exact locations as 
well as the "dynamics" of HAPPI particles. A visual in- 
spection shows that, at high redshift, they are preferentially 
located within the filaments flowing towards halos. At later 
times though, they can be found either in regions of strong 
dynamical activity (e.g., mergers, the outskirts of halos and 
infalling to halos), or at the very centres of relaxed systems. 
The smaller spatial force resolution of the ACDMhappi2 
mentioned earlier can hence be linked to the influence of 
the HAPPI particles at the very centres of halos. 



4 ANALYSIS OF THE MATTER 
DISTRIBUTION 

4.1 Large Scale Structure and Global Properties 

A visual impression of the density field of the particle distri- 
bution at redshift z — for all four models is shown in Fig. 3, 
where the local density at each particle position was de- 
termined by smoothing the distribution onto a regular grid 
(256 3 nodes) and interpolating the density on the mesh back 
to the particle positions. Not surprisingly, the AWDM model 
appears less clumpy and far smoother than the ACDM simu- 
lation. However, there appears to be more small scale struc- 
ture in both the HAPPI runs, or at least the smaller objects 
are more contrasted. We can readily relate this phenomenon 
to the influence of B and larger B values give higher con- 
trasts (remember that the fiducial ACDM model is nothing 
else than a HAPPI model with B=0). But despite the more 
grainy appearance of the HAPPI runs in Fig. 3, the power on 
small scales is reduced compared to the power of the ACDM 
model. This can be verified in Fig. 4, where we plot the dark 
matter power spectrum of density fluctuations for all four 
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Figure 2. The correction term in comparison to the mean-field force acting on a single particle as a function of density. The four panels 
are for all particles in the simulation at redshifts 2 = (upper left), z = 0.5 (upper right), z = 1.0 (lower left), and z = 4.0 (lower 
right). The vertical line indicates the virial overdensity at the respective redshift. The number in each quadrant lists the total number of 
particles in this regime of the plot. 



models at redshifts z = 3 (lower curves) and z = (up- 
per curves). Especially ACDMhappi2 falls behind ACDM 
even though it appears to be marginally more evolved at 
higher redshift. We will though see that these two results, 
the "graininess" of the ACDMhappi2 model and the lack of 
small scale power, do not exclude each other. The absence 
of small scale power on scales below ~ 1-2 h^ 1 Mpc reflects 
the fact that the halos corresponding to those scales (i.e. ha- 
los of mass > 10 n /i -1 M©) are internally less concentrated 
than their ACDM counterparts. 



4-1.1 Beyond the two-point estimators 

Two-point estimators like the power spectrum are sensi- 
tive only to the amplitude of the modes of the fluctuat- 
ing field £>(x) — Qt. The information concerning the rela- 
tive phases is contained in the higher-order correlations. In 



the literature there have been several meaningful quanti- 
ties proposed depending on higher-order correlations with 
a more or less transparent physical interpretation. In this 
work, we have employed the scalar Minkowski Junctionals 
(MFs) (see, e.g., Mecke & Stoyan 2000; Dommguez 2001), 
which allow a quantification of the morphological aspects 
appreciated by visual inspections. Given a density field g(x) 
and a density threshold p, one constructs the isodensity sur- 
face S — {xj£>(x) = p} (with the convention that the re- 
gion q > p is taken as the interior of S). The four MFs 
V v , v = 0, 1, 2, 3, can be defined as surface integrals over 5 
and have the following geometrical meaning (up to a con- 
ventional constant prefactor): 

Vb oc volume enclosed by the isodensity surface S 
Vi oc total area of 5 
V2/V1 oc mean curvature of S, averaged over S 
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Figure 3. Color-coded density field of all four models at rcdshift z = 0. The order is (clockwise starting upper left) ACDM, ACDMhappil, 
ACDMhappi2, and AWDM. 



V3/V1 oc Gaussian curvature of 5, averaged over S 

Vo is in fact proportional to the probability that g(x) > p 
(assuming spatial ergodicity of the realization). The ratio 
V1/V0 is a measure of how compact the volume enclosed 
by 5 is packed. V2/V1 is a measure of the convexity of the 
surface 5 , while V3 is proportional to the Euler characteristic 
or genus of the body defined by S: 

V3 oc number of disconnected objects + number of holes 
— number of tunnels. 



There seems to be a close relation between the threshold 
value at which V3 vanishes and the percolation threshold of 
the volume enclosed by S (Mecke & Wagner 1991; Neher 
2003) — as a matter of fact, the use of percolation analysis 
is not rare in the analysis of cosmological structures (e.g. 
Yess & Shandarin 1996) . As an example of how the MFs are 
to be interpreted, in App. B we discuss the case that £>(x) 
is derived from a realization of a Poissonian distribution of 
points. 

Starting from the positions of the particles given by the 
simulation, a density field £>(x) is generated on a grid by 
smoothing with the window Eq. (12). We generated 20 dif- 
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Figure 4. Dark matter power spectrum P(k) for all models at 
redshifts 2 = 3 (lower curves) and z = (upper curves). 



ferent realizations of the field by displacing randomly the 
grid and we present the MFs averaged over these realiza- 
tions. The MFs are functions of the grid constant and the 
density threshold. We observe that the measurement of the 
MFs in an TV-body simulation is affected by two kinds of 
errors (which are more pronounced for higher orders v of 
the MF K): 

(i) Finite-volume effects: the tails of the probability distri- 
bution of the density field cannot be probed correctly in a 
finite volume, so that the measured MFs are noisy and take 
discrete values (due to the grid) as the threshold approaches 
the extremal values probed in a given simulation. This effect 
can be reduced by increasing the number of realizations. 

(ii) Finite-mass effects: at threshold values corresponding 
to <C 1-10 particles per grid node, the MFs detect that the 
density field is actually derived from a distribution of point 
particles. 

In order to emphasize the differences between the four 
simulations and to facilitate the physical interpretation, we 
plotted as a function of p the relative difference between 
a measure in one model and the corresponding one in the 
reference ACDM model (denoted by 'ref'): 



V v 



rcf 1 (" = 0,3), 



V, 

The plots are shown in Fig. 5 for a grid of 128 3 nodes. 
As a general feature, the trend of the AWDM model with 
respect to the reference ACDM model is always opposite 
to the trend of the two HAPPI runs, with the deviations 
of ACDMhappi2 being larger than those of ACDMhappil. 
Only at very high density thresholds, ACDMhappi2 seems 
to follow an opposite trend to ACDMhappil: this is because 
the maximum density in the ACDMhappi2 run is smaller 
than in the ACDMhappil or ACDM runs. This (physical) 
effect has been already noticed concerning the final force 
resolution of the runs: it is due to the larger effect of the 
additional acceleration C and can be directly linked to the 
location of "HAPPI particles" at redshift z = at the center 
of relaxed structures, see Sees. 3.1 and 3.3. 
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Figure 5. Minkowski functionals vs. density threshold at red- 
shift z = 0. The threshold is given in units of particles per grid 
node. The grid has 128 3 nodes, corresponding to a grid constant 
0.20/i — 1 Mpc, so that the threshold coincides with the ratio g/gt 
for this spatial resolution. 
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The results can be summarized as follows: 

• The dependence Vo(p) quantifies the visual impression 
that low-density regions (g <; 10gb) are clearly less likely 
in the WDM model, while overdense regions (10 < g/gt < 
10 3 ) are more abundant in the HAPPI models. As remarked, 
higher density areas are rarer in ACDMhappi2 and more 
frequent in ACDMhappil compared to ACDM. 

• From the fact that V2, V3 > (not shown in the plot) in 
the whole range of threshold values, one infers that, observed 
at the resolution of 128 3 nodes, the matter distribution con- 
sists mainly of "clusters" (i.e., disconnected objects which 
are convex on average). According to the plot of (V2/Vi)(p), 
these clusters in the HAPPI runs tend to be rounder than 
in the ACDM run, while they tend to be more cigar-shaped 
(filament-like) in the AWDM model. From the plot of Vs(p) 
one deduces that the AWDM model has less clusters at all 
threshold values; the HAPPI runs have more clusters, except 
for ACDMhappi2 at very high densities because this partic- 
ular model does not reach as high densities as the other 
ones. 

• Finally, the plot of (Vi/Vb)(p) shows that matter is 
more compactly packed in the AWDM model, and less com- 
pactly in the HAPPI runs. This is likely due to the different 
cluster abundances just discussed. 

We have repeated the analysis at different spatial resolutions 
(16 3 , 32 3 , 64 3 and 256 3 nodes) for the Minkowski function- 
als. The quantitative differences between models decrease 
as the grid becomes coarser, but the same conclusions hold 
roughly in a qualitative manner. 

Summarizing, the AWDM run has less voids and more 
filamentary-like structures than the reference ACDM model, 
while the HAPPI runs have comparatively more mass con- 
centrated in small, roundish clusters. 



4-1-2 The velocity field 

Motivated by the theoretical discussion, we have also ad- 
dressed the distribution of the (comoving) vorticity, u> = 
V x u, and the (comoving) divergence, 6 = V • u, of the pe- 
culiar velocity field u(x). Using the positions and velocities 
of the particles, we compute the fields w(x) and #(x) on a 
regular grid (see Sec. 3). The cumulative probabilities that 
\uj\ 2 > Cj 1 and that 8 2 > 9 2 are given by the first Minkowski 
functional Vb, which we compute as explained previously. 

Fig. 6 shows P/P Tei — 1, where P is the cumulative 
probability measured in a simulation, and P rcf is the cor- 
responding probability in the reference ACDM simulation. 
The ACDMhappi2 model shows a clear tendency to have a 
larger vorticity and a lower divergence (in absolute value) 
than the ACDM model. The vorticity in the ACDMhappil 
model exhibits a similar tendency, but the differences are 
quantitatively smaller. In order to minimize finite-mass ef- 
fects, we considered only grid nodes for which the value of 
the smoothed density field corresponds to more than 10 par- 
ticle per node. If we restrict the measurement of the proba- 
bility distributions to higher densities, the quantitative dif- 
ferences tend to be smaller. 
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Figure 6. Cumulative probabilities P(\u>\ 2 > u; 2 ), P(B 2 > 6 2 ) 
relative to the probabilities in the ACDM model at redshift z = 
at a resolution of 128 3 grid nodes. The solid symbols represent 
the (relative) cumulative probability of the vorticity, the joined 
open symbols correspond to the divergence. The threshold {Cj 2 , 
8 2 ) is given in units of Hq. 



4.2 Dark Matter Halos 

The analysis in the following Subsection is primarily based 
upon gravitational bound objects which were identified 
using MLAPM's Halo Finder (MHF) (Gill, Knebe & Gibson 
2004). This newly developed halo finder uses the adaptive 
grid structure invoked by the A^-body code MLAPM and re- 
organizes it into a tree. The centres of the grids at the end- 
leaves of a branch of the tree serve as (potential) halo centres 
and all gravitationally bound particles about these centres 
are being collected. For a more elaborate discussion of this 
halo finder we refer the reader to Gill et al. (2004). We only 
like to note at this point that MHF is essentially parameter 
free and naturally finds halos with exactly the same accu- 
racy as the simulation. Besides of the growth of the halo 
mass function and the abundance evolution of gravitation- 
ally bound objects, respectively, we confine our analysis to 
objects identified at redshift z = 0. The investigation of the 
evolutionary history and hierarchical growth of structures 
will be published elsewhere. 



4-2.1 Global Properties 

In Fig. 7 we show the evolution of the cumulative mass func- 
tion n(> M) of dark matter halos, i.e., the number of DM 
halos with mass larger than M. The AWDM model clearly 
has less low-mass objects as already pointed out by other au- 
thors (Knebe et al. 2002, Bode et al. 2001, Colin et al. 2000). 
All three CDM models though perfectly agree at a redshift 
of z = 5, but there is a clear trend for the HAPPI models 
to give rise to more small mass halos, in agreement with the 
conclusion derived in the previous Section. This can be ver- 
ified in Fig. 8, where we plot the number density evolution 
of objects more massive than M > 10 10 /i -1 M (> 20 par- 
ticles). In ACDMhappi2 there are roughly 50% more halos 
above our mass cut than in the reference ACDM run. 

In order to study the clustering properties of these halos 
we estimated the two-point correlation function for low- and 
high-mass objects, respectively, as 
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Figure 7. Cumulative mass function n(> M) of DM halos. 
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Figure 8. Number density evolution of objects more massive 
than M > lO^/i" 1 M . 



?gal(r) = -1 + -^2- ^ 



V n a (r;Ar) 

al 



v(r; Ar 



(15) 



where JV ga i is the total number of objects in the simulation 
volume V, and n a (r; Ar) is the total number of objects in 
a spherical shell of radius r and thickness Ar (and volume 
v(r; Ar)), centered at the a-th object. The result (along with 
Poisson error bars based upon the number of pairs in each 
bin) is shown in Fig. 9. The low-mass objects show very 
similar clustering patterns, with a higher amplitude of £ ga i 
for AWDM and a marginally decreased correlation for the 
ACDMhappi2 model. The situation though is difficult to 
judge at the high-mass end as we have too few pairs per bin 
to make conclusive statements. The errors bars are larger 
than the differences amongst the models and hence have 
been omitted. It seems, however, that the clustering pattern 
of high mass objects is similar in all models and does not 
show differences when including the HAPPI correction term. 
The most striking (and interesting) difference between 
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Figure 9. Two-point correlation function of objects more massive 
(upper panel) and less massive (lower panel) than 10 12 /i -1 Mq. 



standard ACDM and the two HAPPI models, however, 
emerges when we turn to the spin parameter distribution. 
The spin parameter A was calculated using the definition 
given by Bullock et al. (2001a), 



A = 



V2AL 



(16) 



vir <Aar ' vir 



where L is the angular momentum of the halo with respect 
to its center of mass, r v i r is the virial radius of the halo, M v i r 
is the virial mass (mass enclosed within the virial radius), 
and tVir = J GM v ir/r v i r is the circular velocity at the virial 
radius. The virial radius and mass are determined by the 
condition 



47T 3 

M vil = — r vir £ v 



(17) 



where q v [ t = A v i r Qb(z = 0) is a fiducial density with A v i r ~ 
340 (at redshift z = 0) based on the dissipationless spherical 
top-hat collapse model for the cosmological parameters of 
the ACDM model. The probability distribution, P(A), of 
the spin parameter was fitted to a log-normal distribution 
(e.g. Frenk et al. 1988; Warren et al. 1992; Cole & Lacey 
1996; Mailer, Dekel & Somerville 2002; Gardner 2001), 
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Figure 10. Upper plot: Measured spin parameter distribution 
for all four models at redshift z = 0, and the corresponding fit 
to Eq. (18). Lower plot: Measured cumulative distribution of the 
spin parameter. 



P(X) = 



exp 



ln 2 (A/A ) 
2ol 



(18) 



The results are presented in Fig. 10 and in Table 2, show- 
ing that a larger value of B entails a larger spin parameter. 
In the lower panel of Fig. 10 - where we plot the cumula- 
tive distribution of the spin parameter - we clearly see that 
for a given A the probability to find halos with a smaller A 
is greatly lowered in ACDMhappi2 — or in other words, it 
is more likely to find halos with larger spin parameters in 
ACDMhappi2. As we will see later on (cf. Fig. 13 in Sec- 
tion 4.2.2), there is a mass dependence of this result: lower 
mass halos tend to dominate the signal seen in Fig. 10. When 
plotting the mass-weighted spin parameter distribution (not 
shown) the peaks for all four models approach each other. 
However, even then the HAPPI runs still show a distinct tail 
out to larger A values. 

The HAPPI correction term also affects the concentra- 
tion of dark matter halos. We define the concentration C1/5 
as the ratio of the virial radius r v i r and the radius of the 



Table 2. Parameters derived from fitting the spin parameter dis- 
tributions P(X) to Eq. (18). 



model 



Ao 



ACDM 0.0287 0.4485 

ACDMhappil 0.0333 0.4611 

ACDMhappi2 0.0596 0.4914 

AWDM 0.0259 0.4001 




ACDM 

ACDM happi 1 
ACDM happi 2 
AWDM 



0.001 -4 



Figure 11. Cumulative distribution of the concentration param- 
eter c 1/5 . 



sphere that contains 1/5 of the virial mass (i.e. ri/5 is de- 
fined via M(< n/ 5 ) = Mvir/5): 



C1/5 = 



n/5 



(19) 



In view of the definition (17), it follows that the average 
density within a radius rys is given by (l/5)ci/ 5 £? v ir- Fig. 11 
plots the cumulative probability distribution of the concen- 
tration of halos. We observe an obvious trend for an over- 
abundance of low-concentration halos in the AWDM and 
ACDMhappi2 models. However, the opposite actually holds 
for ACDMhappil, where there appear to be of order 10% 
more concentrated halos. The relative lack of power on scales 
<; l/i" 1 Mpc noted in Fig. 4 for WDM and ACDMhappi2 is 
related to the relatively lower concentration (and increased 
smoothness for WDM) of the halos observed in these models 
(one must bear in mind that the halos have a virial radius 
< 800ft" 1 kpc, see Fig. 16). 

Although not shown, we confirm that the scaling of the 
concentration C1/5, defined by Eq. (19), with mass M v i r 
follows the same relation as the one proposed by Bul- 
lock et al. (2001b) for the "NFW concentration" of the halo, 
cnfw = r v i T /r a (see Eq. (21) for the definition of r s ). 



4-2.2 Cross Correlations 

While the last Section dealt with the distribution of halo 
properties, we now compare these properties across the mod- 
els, i.e. how do these properties change in a given halo when 
moving from one model to another? 
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In order to find corresponding halos across the four 
models, we compare their individual particle content. We 
start with a halo in the ACDM model, whose particles are 
tagged, and locate the corresponding halo in the other three 
models as the one that shares the largest number of tagged 
particles. In the Figs. (12-14) we always plot the value of 
the property under investigation in the ACDM model, di- 
vided by the value of said property in the other model for 
all "cross-identified" halos. These "scatter-plots" are always 
accompanied by histograms, where we average the ratios 
presented in the respective figure in nine bins across the ac- 
tual mass range. The percentages of counterparts amount to 
practically 100% for the HAPPI models while there are only 
40% cross-identified halos in the AWDM model. This num- 
ber though increases to about 90% when we consider halos 
containing more than 200 particles (i.e., M> lO 11 /! -1 Mg), 
and hence we set this as a lower limit in the cross-correlation 
plots. 

We start with the most obvious halo attribute, namely 
the halo mass itself. Fig. 12 shows that there is a very tight 
correlation for the masses of individual halos, especially for 
ACDM and ACDMhappil. The scatter about the 1:1 re- 
lation (the flat line of value 1) marginally increases from 
~11% (averaged la-value) for ACDMhappil to ~18% for 
ACDMhappi2. At the high-mass end, the ACDMhappi2 ha- 
los tend to have a slightly larger mass. The most pronounced 
differences can be found for AWDM though: at the low-mass 
end the halos in AWDM appear to be less massive than their 
ACDM counter parts. 

In the previous Section we showed that the HAPPI cor- 
rection term lead to an increase in angular momentum by 
investigating the spin parameter distributions. But can we 
be sure that the observed rise in A as defined by Eq. (16) is 
not related to a possible decrease of virial radius r v i r and/or 
virial velocity v V ir? To clarify this uncertainty we show in 
Fig. 13 the cross correlation of total specific angular mo- 
mentum, 
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(20) 



In view of the already mentioned minimal scatter in the mass 
of cross-identified halos between the ACDM model and the 
HAPPI models, Fig. 13 confirms the previous result that one 
effect of the term C in Eq. (7b) is to inject angular momen- 
tum to halos. We also note that there is a mass dependence 
in this trend: the differences in angular momentum are on 
average larger for lower mass objects, this being particu- 
larly noticeable for the ACDMhappi2 and AWDM models; 
this "break" roughly happens at around 10 12 /i _1 M©. 

We close this Section with an investigation of the cross- 
correlation of the concentration parameter C1/5 in Fig. 14. 
The results are consistent with Fig. 11: ACDMhappil halos 
are more concentrated than their ACDM counterparts, while 
the excess of low-concentration halos for ACDMhappi2 is 
due to high-mass halos. The mass trend already noted in 
Fig. 13 can also be acknowledged in this figure. 

Finally, we mention as a general property that the dis- 
persion in the scatter plots increases as the halo mass dimin- 
ishes. We attribute most of this scatter to differences in the 
halo's formation history, but a detailed study as function of 
redshift is required which we will postpone to a later paper. 



Figure 12. Ratio of halo masses for cross-identified objetcs. The 
virial mass M v ; r is given in units of h _1 M©. The histograms 
represent the mean ratio in the respective bin. 



Moreover, numerical effects could also contribute to some 
extent. 



4-2.3 A Closer View of Individual Halos 

A visual representation of the two most massive halos in 
all four models is given in Fig. 15. This figure nicely demon- 
strates the result regarding the lower concentrations in (high 
mass) ACDMhappi2 halos: the second most massive halo 
does not even show a distinct centre in ACDMhappi2 and 
appears more "puffy" than in any of the other models. 

The question now arises whether the density profiles of 
the HAPPI models can still be fitted by the universal density 
profile advocated by Navarro, Frenk & White (NFW, 1997) 



p{r) 



PsT s 



r(r s + r) 2 



(21) 
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Figure 13. Ratio of specific angular momenta. 



Figure 14. Ratio of halo concentrations C1/5. 



Fig. 16 now shows p(r) and the corresponding best fits to 
NFW profiles for a selection of halos covering the mass range 
from the most massive one (upper left) to rather light halos 
(lower right) containing a mere 300 particles. For (nearly) 
all ACDMhappi2 halos we observe a relative flattening in 
the central regions. 

In order to gauge the quality of the fits (for the 16 
presented sample profiles) we calculate the \ value defined 



2 

X = 



N hins 



iv bina 



1=1 



PNFw(?"i) 



PNFW^i 



(22) 



where pi are the binned density profiles derived from the 
simulation data and Pnfw the best-fit NFW profiles. This 
analysis then indicates that all four models are equally well 
fit by Eq. (21) with x varying in the range \ 2 ~ 0.02-0.05 
depending on the weighing scheme applied for each indi- 
vidual bin. This entails that the dark matter halos of the 



HAPPI runs still exhibit the rather infamous "cusp" at the 
center. 

We remind the reader again that the force resolution 
throughout the runs varies. Whereas ACDM, ACDMhappil, 
and A WDM reached 2.5/i _1 kpc resolution, ACDMhappi2 
reliably resolves structures only on scales larger than 
10/i -1 kpc. Moreover, the resolution can also change from 
halo to halo due to the adaptive mesh nature of both the 
halo finder and the iV-body code: not all halo centres lie on 
the finest grid level reached in the simulation. However, we 
plot profiles starting from the distance r m i n that corresponds 
to a sphere containing at least 10 particles (and hence r m i n 
can be actually smaller than the nominal resolution of the 
simulation) . 

For the same set of halos we present in Fig. 17 the 
rotation curves out to half the respective virial radius. The 
rotational velocity v c i IC (r) is defined as 
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Figure 15. Visual representation of the two most massive halos 
at redshift z = 0. The left column shows the most massive halo 
while the right column is the second most massive one. From 
top to bottom one has ACDM, ACDMhappil, ACDMhappi2, and 
AWDM. 



«circ(V) = 



2 GM{< r) 



(23) 



There are a number of interesting observations to discuss 
now. We find that in nearly every halo the ACDMhappil 
rotation curve rises to higher values than any of the other 
models. While the maximum is still at comparable distances 
in ACDM and ACDMhappil, the latter shows a steeper in- 
ner increase and a subsequent faster decline to nearly the 
same level in the "outer" parts. Moreover, the ACDMhappil 
rotation curves are always slightly above the corresponding 
ACDM curves. Quite the opposite is true for ACDMhappi2. 
Here we find that in most of the cases the circular rota- 



tion values at a given radius are substantially smaller than 
in ACDM. However, this difference becomes less prominent 
in lower mass systems and the flat part of the rotation 
curve nearly reaches the same level as ACDM. This discrep- 
ancy in the trends between ACDMhappil and ACDMhappi2 
for high-mass halos is likely related to the also opposite 
trends concerning the concentration, Figs. 11 and 14, and 
the small-scale power, Fig. 4. 



5 DISCUSSION & CONCLUSIONS 

We have presented a series of cosmological A^-body simula- 
tions which made use of the hydrodynamic approach to the 
evolution of structures (Dommguez 2000) . This approach is 
novel in that it deals with the mass density and velocity 
fields with explicit account of the coarse-grained nature in- 
trinsic to any approach of solving, for instance, Poisson's 
equation via Monte Carlo sampling of phase-space. This N- 
body approach unavoidably introduces finite resolution ef- 
fects and there have been systematic studies of the conse- 
quences in the context of cosmological structure formation 
(Kuhlman, Melott & Shandarin 1996; Splinter et al. 1998; 
Moore et al. 1998; Knebe et al. 2000; Power et al. 2003). N- 
body simulations invariably neglect the dynamical effect of 
sub-resolution degrees of freedom altogether. For the first 
time we have run simulations including a physical model 
of the coupling to the neglected scales. A^-body codes are 
usually viewed as integrators of the Vlasov-Poisson system 
of equations. However, we have argued how grid-based N- 
body codes such as MLAPM can be reinterpreted to integrate 
hydrodynamic-like equations for the mass density and ve- 
locity fields. 

The additional, correction term introduced in the hy- 
drodynamic approach is proportional to a "coupling con- 
stant" B which depends on the smoothing window used to 
calculate the coarse-grained fields. It is found to be B — 1/4 
for the triangular-shaped-cloud window used throughout 
the A^-body code MLAPM. In order to get a better understand- 
ing of the effects of the correction term onto the evolution 
of cosmic structures we also performed a simulation with 
a higher value B = 1 — this later model is not physically 
motivated but rather serves as an "academic toy model" for 
comparison. The standard ACDM simulation can be under- 
stood as another HAPPI run with the value B — 0. In order 
to allow for a better comparison with other feasible alterna- 
tives to the concordance ACDM model as well as to better 
gauge the influence of the correction term, we also simulated 
the evolution of the same structures in a AWDM universe. 

In this work we concentrated on the comparison of the 
four simulations at redshift z — 0. We analyzed the resulting 
structures in two complementary manners: global properties 
of the mass density and velocity fields, on the one hand, and 
specific properties of DM halos, on the other hand. We find 
appreciable differences between the B ^= runs and the ref- 
erence (B — 0) ACDM run, even though the force due to 
the correction terms are for most particles one or even two 
orders of magnitude smaller than the total force (cf. Fig. 2). 
Most remarkably, the correction term favors the prolifera- 
tion of low-mass halos, giving the mass distribution a more 
"grainy" aspect, as well as the gain of angular momentum 
specially by low-mass halos, which also shows up in a ve- 
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Figure 16. Density profiles for halos at redshift z = 0. The number printed into each halo panel is the mass of the halo in units of 
h~ 1 Mq. The vertical lines to the right indicate the respective virial radius while the two vertical lines to the left indicate the spatial 
resolution of the ACDM, ACDMhappil, A WDM (dashed) and ACDMhappi2 (dotted) run, respectively. 



locity field with a larger vorticity. These effects are quanti- 
tatively larger as the value of B increases; for B = 1/4 the 
differences lie at the (10 — 20)% level (and even higher for 
the specific angular momentum at low masses). A feature 
in which the B = 1/4 and B — 1 runs exhibit an opposite 
trend with respect to the B = run is the concentration of 
high-mass halos: the B = 1 run results in an overabundance 
of high-mass halos with a lower concentration This is par- 
alleled by a smaller circular velocity of these halos, and by 
the relative lack of power in the spectrum of density fluc- 
tuations at sufficiently small scales, so that the maximum 
density reached in the B = 1 run is much smaller than in 
the other runs. The B — 1/4 run, however, shows precisely 
the opposite tendency with respect to the reference run. One 
can conjecture that this discrepancy between the B = 1/4 
and 5=1 runs lies in a difference in the rate of shear and 
vorticity generation and of kinetic energy drainage by the 
correction term. A comparative study of the structures at 
different redshifts is required in order to obtain more precise 
conclusions about this issue. 

The relatively small quantitative differences between 
the B = 1/4 and the 5 = runs evidenced in the prop- 



erties that we have measured suggest that the B = 1/4 
correction term could be considered a small perturbation to 
the B = evolution. By contrast, the results of the B = 1 
run indicate that the correction term should not be treated 
as a perturbation in this case. 

Our results agree with the theoretical expectation for 
the qualitative behavior of the correction term, which mod- 
els the effect of small-scale tidal torques and shear stresses 
(Dommguez 2000, 2002; Buchert & Dommguez 2005). We 
observed that the correction term is dominant preferentially 
in walls at high redshifts, and later on in filaments, regions of 
mass accretion onto halos, halo centers as well as in regions 
of particular dynamical activity (e.g., mergers), that is, re- 
gions of large gradients in the fields, in concordance with the 
form of the correction term (8). The term is expected to act 
as a drain of kinetic energy in collapsing regions: this can 
explain the formation of small clusters of particles, which 
would otherwise fly by each other — instead, they can be 
gravitationally confined by a potential well that is lower than 
in the B = model. This could explain the low mass halo 
proliferation for B = 1/4 or B = 1, as well as the observed 
tendency of halos to attain a slightly more concentrated con- 
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Figure 17. Rotation curves for the same halos as in Fig. 16. 



figuration in the case of B = 1/4, when the correction term 
can be considered a small perturbation. If B = 1, on the 
other hand, the loss of kinetic energy is apparently so im- 
portant that, in some cases of halos in regions of high dy- 
namical activity, dynamical relaxation and coalescence are 
slowed down noticeably, leading to a multiple-core struc- 
ture. These not completely relaxed halos would then have 
a lower concentration and a lower mass than their LCDM 
counterparts, similarly to the simulation results. 

We further confirmed explicitly that the correction term 
acts as a source of vorticity. This relates directly to the gain 
in angular momentum of halos, which tends to increase with 
growing value of B, especially at the low-mass end of the 
halo distribution. 

Finally, we remark that our findings agree qualitatively 
with conclusions following from a comparative study of iden- 
tical initial conditions evolved at different resolutions. We 
have run a series of test simulations where we either switched 
on the HAPPI correction term or increased the actual force 
resolution; both methods lead to comparable results that 
are in qualitative agreement with the conclusions presented 
here. For a more quantitative analysis we though refer the 
reader to a future paper in preparation where we will in- 
vestigate the relationship between HAPPI simulations and 



higher-resolution ones in more detail. The proliferation of 
small halos with increasing resolution has also been reported 
by other authors in and around (massive) halos (Klypin et 
al. 1999, Moore et al. 1998) as well as in voids (Gottlober 
et al. 2003). Concerning the generation of angular momen- 
tum though, the relevance for the formation of realistic disk 
galaxies has yet to be determined but there are clear indica- 
tions that this task requires good mass and force resolution 
(Governato et al. 2004). In conclusion, the HAPPI imple- 
mentation seems indeed to be qualitatively consistent with 
what one expects from higher resolution simulations and 
hence may provide a framework for a better understanding 
of resolution effects in iV-body simulations. However, further 
work is required to substantiate this possibility. 
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APPENDIX A: CONSERVATION OF ENERGY, 
MOMENTUM AND ANGULAR MOMENTUM 

The coupling to the small-scales modelled by C in Eqs. (7) 
injects energy into (or drains energy from) the resolved spa- 
tial scales. Given the simulation box Vbox with periodic 
boundary conditions, we can define the total peculiar kinetic 
energy and the total peculiar mean-field potential energy as 
follows: 



K 



r 



dx gu , 



(Al) 



1 2 
« a 

2g b 



dxdy [g(x) - Qb][e(y) - Qb]S(x - y), 



V bo 



where the time-independent, symmetric kernel S(x) is the 
solution of the problem 

V 2 5 = 47rG(a 3 £l6 ) ( 5 Dirac (x) (x G V box ). (A2) 
The mean-field gravitational acceleration is given by 



w mf (x) = - -^V x 



dy[0(y)-0b]S(x-y). (A3) 



Then, it is easy to show from Eqs. (7) that the total peculiar 
energy H = K + U satisfies the evolution equation 



~dt 



= -H(2K + U' 



Jv bo 



dx gu ■ C. 



(A4) 



This is a generalization of the Layzer-Irvine equation. Due 
to the correction term, the condition of "mean-field virial- 
ization", 2K-\-U lni = 0, does not imply a time-independent 
H. The quantity 

1= aH + 



J daK 



(A5) 



is conserved by the original Layzer-Irvine equation but is 
not constant according to the generalized equation (A4). 

Concerning momentum and angular momentum, 
Eqs. (7) do not violate global conservation. Let V(t) denote 
a time-dependent volume defined by the condition that the 
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mass enclosed is constant, i.e., a Lagrangian domain. The 
peculiar momentum of the domain, 



Vv 



dx gu, 



V(t) 



verifies the evolution equation 
dVv 



dt 



= -HVv 



dx g(\ 



(A6) 



(A7) 



V{t) 



The correction term can be written as the divergence of a 
tensor (Buchert & Domfnguez 2005) 



gCi = Bl/dj 



g(diwf f ) + 2%GaQ 2 8ij g(d k Ui){d k Uj) 



so that its contribution in Eq. (A7) is a surface integral over 
the border of V(t). In particular, when V(t) = Vbox, this 
surface integral vanishes by periodic boundary conditions 
and, since the contribution by w mf also vanishes in this case, 
Eq. (A7) states that aPv box is a constant of motion. 

In the same manner, one defines the angular momentum 
of the domain V(t) with respect to its center of mass X cm (£), 



dx (x - X cm ) x gu. 



V{t) 



The evolution equation for this quantity is 



dCv 
dt 



dx (x — X c 



X £>(w mf ' + C). 



(A8) 



(A9) 



V{t) 



The contribution by C can be written again as a surface 
integral over the border of V(t). Thus, when V(t) = Vbox, 
Eq. (A9) predicts that Cv hox is also a constant of motion. 

As discussed in Sec. 2, the correction C is a source of 
vorticity in the otherwise curl-free flow of the " dust model" . 
Eq. (A9) shows that the correction also affects the evolution 
of the angular momentum of a domain. In this case, how- 
ever, the effect may not be so noticeable, since already at 
the level of the "dust model" there are tidal torques by the 
mean-field gravity w mf . Moreover, since the contribution by 
C is a surface integral, it may be expected to be less relevant 
for a larger domain V(t). Actually, we can rewrite the defini- 
tion (A8) by inserting the identity 2(x— X cm ) = V|x— X cm | 2 
as 



-a 4 <£> dS x £>u|x — X cm | 2 

BV 



(A10) 



ia 4 y dx jx - X cm | 2 [guj - (V g) X u], 

and we see that vorticity is but one contribution to the an- 
gular momentum of a Lagrangian domain. 
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Figure Bl. Minkowski functionals vs. density threshold of the 
realization of a Poissonian distribution of points. 

• The volume Vo(p) decreases monotonically as the 
threshold is increased and the high-density regions (g > p) 
shrink. 

• The area Vi(p) first increases as the low-density regions 
(g < p) expand and, after reaching a maximum, it decreases 
as the high-density regions (g > p) shrink. 

• The average mean curvature V2(p) increases 
monotonously from a negative value (5 is concave towards 
the shrinking high-density region) to a positive value (5 is 
convex towards the shrinking high-density region) . 

• Finally, the genus Vs(p) is positive when 5 looks 
bubble-like: there are many unconnected expanding low- 
density regions ( "holes" ) at small p, and many unconnected 
shrinking high-density regions ("clusters") at large p. Vss 
is negative when S is predominantly saddle-shaped: one 
observes many intertwined high- and low-density regions 
("tunnels") at intermediate p. 



APPENDIX B: MINKOWSKI FUNCTIONALS 
OF A POISSON DISTRIBUTION 

As an illustration of the dependence of the MFs on the 
threshold, Fig. Bl shows the MFs of a realization of a Pois- 
son distribution of points: 128 3 particles were distributed 
randomly in a cubical box, and the density field g(x) was 
obtained by smoothing with the window Eq. (12) in a cu- 
bic grid of 16 3 nodes. The plots are symmetric about the 
mean value of the density (512 particles per node) and span 
a width w rms density (= v512) along the threshold axis. 
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